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£ ABSTRACT 

i We present a series of N-body experiments which confirm the reality of the previrializa- 

tion effect. We also use weakly nonlinear perturbative approach to study the phenomenon. 

^ These two approaches agree when the rms density contrast, a, is small; more surprisingly, 
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they remain in agreement when a ?s 1. When the slope of the initial power spectrum is 
O n > — 1, nonlinear tidal interactions slow down the growth of density fluctuations and the 
On magnitude of the suppression increases when n (i.e. the relative amount of small scale 
^ power) is increased. For n < — 1 we see an opposite effect: the fluctuations grow more 
^ rapidly than in linear theory. The transition occurs at n = — 1 when the weakly nonlinear 
£j correction to a is close to zero and the growth rate is close to linear. Our results resolve re- 
c3 cent controversy between two N-body studies of previrialization. Peebles (1990) assumed 
n = and found strong evidence in support of previrialization, while Evrard & Crone 
(1992), who assumed n = —1, reached opposite conclusions. As we show here, the initial 
conditions with n = — 1 are rather special because the nonlinear effects nearly cancel out 
for that particular spectrum. In addition to our calculations for scale-free initial spectra, 
we show results for a more realistic spectrum of Peacock & Dodds (1994). Its slope near 
the scale usually adopted for normalization is close to —1, so a is close to linear. Our 
results retroactively justify linear normalization at 8/ir 1 Mpc, while also demonstrating 
the danger and limitations of this practice. 

Subject headings: cosmology: theory - galaxies: clustering - galaxies: formation - 
large-scale structure of universe 

To appear in ApJ 



1 



1 Introduction 



It is generally believed that the structures observed in the universe today grew by gravita- 
tional instability out of much smaller inhomogeneities present at the epoch of the decou- 
pling of matter and radiation. The details of this general picture, however, are still under 
discussion. One of the sources of disagreement is the size of the so called previrialization 
effect. The term "previrialization" was introduced by Davis & Peebles (1977), although 
the general idea was expressed earlier (Peebles & Groth 1976). 

According to the previrialization conjecture, small scale fluctuations present within col- 
lapsing mass concentrations tend to slow down their growth. Initial asphericities and tidal 
interactions between neighboring density fluctuations induce significant nonradial motions, 
which oppose the collapse. Virialized clumps form later than predicted by frequently used 
approximations, such as linear perturbation theory or the spherical collapse solution (an- 
other way to state the problem is to say that to obtain a given final density contrast, in 
case of a realistic mass distribution, one needs a significantly larger initial contrast than the 
one for an isolated spherical fluctuation). Some N-body experiments (Villumsen & Davis 
1986; Peebles 1990) appear to reproduce these effects, while others do not. An example of 
the latter is the paper by Evrard & Crone (1992). 

Using his least action method of evolving particle orbits back in time, Peebles (1990) 
compared histories of isolated clumps with that of an ensemble of mass concentrations. The 
required initial density contrast in the former case was consistent with the spherical model. 
The initial contrast required for the more realistic mass distribution was significantly larger. 
Although the spherical model predicts that the growth of the density contrast is faster than 
the linear approximation, the growth observed in the numerical action method is slower. 

Evrard & Crone (1992) posed the question, whether small-scale structure affects the 
clustering on larger scales. Their N-body simulations supported the conclusion that the 
abundance of rich clusters finally formed was insensitive to the amount of small-scale power 
present in the initial conditions. One of our purposes here is to address this controversy. 

Interactions between fluctuations at different scales can be also studied using weakly 
nonlinear perturbation theory. A formalism, appropriate for this purpose, was proposed 
over a decade ago (Peebles 1980; Juszkiewicz 1981; Vishniac 1983; Fry 1984; Juszkiewicz, 
Sonoda & Barrow 1984). This formalism uses Newtonian hydrodynamics and an important 
assumption that the linear term in the perturbative expansion for the density contrast 
is a Gaussian random field. Unfortunately, N-body simulations available in the early 
eighties were too crude to study the range of validity of perturbative methods. This made 
further progress difficult and the field became dormant for several years. With time, the 
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dynamical range of the simulations improved. Moreover, in the last decade the depth of 
galaxy surveys has increased by an order of magnitude, allowing accurate measurements 
of spatial correlations of galaxies at separations large enough to make weakly nonlinear 
theory useful in their analysis. 

As a result of these developments, the study of the evolution of the mass distribution 
in the expanding universe has become a rich topic for analytic perturbative methods. The 
renewed enthusiasm for analytic calculations produced many new results (for a review, 
see e.g. Bouchet & Juszkiewicz 1993; Sahni & Coles 1995; Juszkiewicz & Bouchet 1996 
and references therein). Those directly relevant for our purposes here are contained in the 
papers of Suto & Sasaki (1991) and Makino, Sasaki & Suto (1992). Studying the nonlinear 
evolution of the power spectrum, these authors succeeded in reducing the dimensionality 
of the mode coupling integrals, which appear in the formalism of Vishniac (1983) and 
Juszkiewicz et al. (1984). As a result, they were able to obtain closed-form expressions 
for spectral corrections for scale-free initial conditions. Using the same formalism, Jain &: 
Bertschinger (1994) investigated the transfer of power in the CDM spectrum. They com- 
pared their results to N-body simulations, finding a good agreement in a wide dynamical 
range. 

In this paper, we address the previrialization controversy using weakly nonlinear pertur- 
bation theory. As a measure of the previrialization effect we calculate the weakly nonlinear 
corrections to the second moment (the variance) of the cosmic density field. As it was al- 
ready pointed out by Juszkiewicz et al. (1984), and confirmed by Suto & Sasaki (1991), 
Makino et al. (1992), and Jain & Bertschinger (1994), the strength and the wavenumber 
dependence of nonlinear corrections to the power spectrum are determined by the shape of 
the initial spectrum. Hence, we may expect that the strength of the previrialization effect 
may also vary with initial conditions. As we will show here, this is indeed the case and 
the differences in conclusions, reached by Peebles (1990) and Evrard & Crone (1992) may 
actually reflect the differences in the slopes of the power spectra in their simulations. 

This paper is organized as follows. In section 2, we summarize the formalism of weakly 
nonlinear perturbation theory and calculate the second order contributions to the variance. 
We use Gaussian initial conditions and perform the calculations for scale-free spectra and 
a more realistic spectrum of Peacock & Dodds (1994), based on data from several galaxy 
surveys. The results and discussion are given in Section 3, where we check our perturbative 
calculations against a set of N-body simulations. We also compare perturbation theory 
with calculations, based on a phenomenological model of nonlinear clustering, similar to 
that of Hamilton et al. (1991). 
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2 Perturbation theory 



We assume a model universe with vanishing cosmological constant and arbitrary density 
parameter Q. The content of the universe is supposed to behave as a pressureless fluid 
which undergoes gravitational evolution described by the usual Newtonian equations. The 
cosmic density field is characterized by the density contrast S = S(x,t) = Sp/pb, where x 
is the Eulerian comoving coordinate, t is the cosmological time, pb denotes the background 
density, and Sp is the difference between local mass density and pb- 

The perturbative expansion of the density contrast around the background solution 
6 = is 

6 = S 1 + <5 2 + <$ 3 + • • • (1) 

where S n = 0(6±), Si being the linear theory solution. The n-th order solutions are 
obtained from equations describing the Newtonian evolution using the solutions of the 
(n — l)-th order of density and velocity fields as source terms (see Fry 1984; Goroff et al. 
1986). 

In the Einstein-de Sitter universe the scale factor scales with time as a a £ 2 / 3 and the 
background density pb = 3a 2 /87rGa 2 = l/QnGt 2 . The time dependence of the n-th order 
follows 

6 n (x,t) = [D(t)] n 6 n (x) (2) 

where D(t) oc a(t) and we consider only the mode growing in time. For an arbitrary 
cosmological model, however, the time dependences of different orders should be considered 
independently. Fortunately for a wide range of Q the solutions for the density contrast are 
very weakly dependent on Q and the Einstein-de Sitter case provides a good approximation 
(Bouchet et al. 1992, 1995; Catelan et al. 1995). 

All of the following calculations are much simpler if they are performed in Fourier space. 
For the first order of the density contrast field we have 

<5i(k, t) = D{t) J d 3 x ^(x)e- ik - x (3) 

and the inverse Fourier transform is 

«i(x, t) = D(t)(2n)- 3 J d 3 k ^(k)e ik - x . (4) 

For the following calculation only second and third order solutions for the density 
contrast are needed, and we give them here in the Fourier representation (e.g. Goroff et 
al. 1986). 

S 2 (^t) = -^- J d 3 p J d 3 g^(p + q-k)5 1 (p)^ 1 (q)P 2 (s) (p,q) (5) 
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« 3 (M) = J d 3 p J d 3 q J d 3 m ^(p + q + m-k)^ 1 (p)5 1 (q)^ 1 (m)P 3 (s) (p,q,m). (6) 



(2*) 

The symbol 6u represents the Dirac delta. The symmetrized kernels are of the form 

A (s) (p>q) = :j^(p + q,p,q) 



(7) 



^i S) (P,q ? m ) = A [ #(p + q + m,p)J(q + m,q,m) + 

+ i?(p + q + m, q + m) L(q + m, q, m)] + 
+B F(p + q + m, p, q + m)L(q + m, q, m) + 
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where A = 1/108 and B = 1/189. In the expression above the notation follows that of 
Makino et al. (1992) i.e. 



#(p,q) 

^(p + q, p, q) 
^(p + q,p,q) 
£(p + q,p,q) 



p q 



1 |p + q| 2 pq 

2 p 2 q 2 

4 (E^! + 7^p.q + 10 



8 (E^ + 7 p!±j! p . q + 6 . 



p 2 q 2 



p 2 q 2 



(9) 
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(11) 
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The smoothing of the fields is introduced by the convolution of the density and the 
filtering function 

6 R (x,t) = Jd 3 y6(y,t)w(\x-y\,R) . (13) 

We use lower-case w for spatial filters and capital W for their Fourier representation, given 
by 

W(k) = [ d 3 xw(x)e- ik -* . (14) 



(15) 



Our window functions are normalized so, that 

J d 3 xw(x) = W(0) = 1 . 



We perform our calculations for two types of windows - a Gaussian and a top-hat with 
Fourier representations respectively 



W G (kR) = e- k2R2 ' 2 



(16) 
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and 



W TH (kR) = 3^(kR)-*Jz(kR) (17) 



where J3 is a Bessel function 
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= £ - C ° S(I) ) ■ (18) 

We assume a Gaussian distribution for 61 in equation (1) and define 

cr 2 = <«5 2 ) = D\t) J ^ P(*) IT 2 ^) (19) 

as the linear variance of the density field. We assume that for cr < 1, the first few terms 
in the perturbative expansion provide a good approximation of the exact solution (this 
assumption is supported by results of N-body simulations, discussed below). Since Si is 
assumed to be a Gaussian random field, all its statistical properties as well as those of the 
higher order terms in the perturbative series (1) are determined by the power spectrum 
P(k), defined as 

(^i(k)^(q)) = (2nfS D (k + q) P(k) . (20) 
We begin by considering spectra with a power-law form 

P(k) = Ck n , -3 < n < 1 , (21) 

where C is a normalization constant. For such fields, smoothed with a Gaussian filter, the 
linear order contribution to the variance (equation 19) is 

TV n+3 \ 

4 = CDHt)^^ • (22) 

For a top-hat smoothing we get 

a 2 TH = CD\t) 87r3 /2^ + 3r(i _ f)r( V) ' ^ 

The perturbative series (1) can be used to expand the second moment of the density 
field in powers of cr, 

(^ 2 ) = (Si) + (<5 2 ) + 2(M 3 > + • • • = a 2 + (/ 22 + / 13 y 4 + 0(a 6 ). (24) 

The first term is the linear variance given by equation (19). The two following terms are the 
lowest order nonlinear corrections. These terms are of order cr 4 because Si is a Gaussian 
random field (all odd-order terms in the expansion, including those of order cr 3 , must 
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vanish). We have introduced the quantities J 2 2 and ii 3 because it is convenient to calculate 
the first nonlinear corrections in a normalized form, i.e. divided by cr 4 . Indeed, provided 
the integrals are convergent, for the scale-free power-law spectra, these ratios should be 
dimensionless numbers, independent of scale, as in the case of the higher moments (e.g. 
the skewness and the kurtosis) of the fields. By using the second and third order solutions 
(5) and (6) we obtain 



D\t) 
22 98(2vr)V 



J d 3 p J d 3 qP(p)P(q)W 2 (\p + q\R) J 2 (p + q,p,q) (25) 



ha = j d?p j d?q P{p)P{q)W\qR) 



x { A[ H(q,-p) J(p + q,p,q) 

+#(q,p + q) L(p + q,p,q)] (26) 

p -> q 



B F(q, -p, p + q) L(p + q, p, q)} + 



q -> p 



where the last term in brackets means that similar expression with p and q interchanged 
should be added if the symmetry in p and q is to be maintained. After integration the 
expressions become 

D\t) 



iij = r dk p w2 ( kR ) p ^) > ( 2? ) 



where 



*4 - 158 + 100x 2 - 42x 4 + ^-(x 2 - l) 3 (7x 2 + 2) In 1 + ' 



(29) 



The window function in equation (27) describes spatial smoothing, applied to the density 
field; each term in the perturbative expansion (1) was convolved with a Fourier transform 
oiW(kR). 

The input from I22 is always positive and from /13 negative, therefore we may inter- 
pret their values at some scale respectively as the additional power coming in from other 
wavelengths and as the power that is lost and taken over by other wavelengths. 
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Both J 2 2 and I i3 diverge individually in the limit of k — > and kx = q — > if n < — 1. 
Fortunately, in the entire interval —3 < n < 1, all diverging terms cancel when I22 and /13 
are added together, so that their sum, 

I 2 = /22 + /13 (30) 

remains finite. These "miraculous cancellations" are probably associated with the Galilean 
invariance of the Newtonian equations of motion (Scoccimarro & Frieman 1995). In the 
opposite limit of k — > 00, q — > 00, divergencies occur for J 2 2 if n > \ and for 7 13 if n > — 1. 
Therefore if we want to calculate the sum of both terms for n > —1 we need to introduce 
a cutoff in the initial power spectrum at large wave-numbers 

I for > & c . 

An equivalent way to introduce the cutoff is to smooth the initial power spectrum with a 
Gaussian 

P(k) = Ck n e- k2r2 . (32) 

We introduce shortwave cutoffs in the power spectrum to avoid divergencies in pertur- 
bation theory for scale-free spectra. Diverging terms in the perturbative expansion do not 
necessarily imply that the fully nonlinear solution for a(R,t) is divergent; strong nonlin- 
earities may "naturally" introduce a nonperturbative cutoff. Moreover, spectra without 
a cutoff represent a rather far-fetched idealization anyway. In any realistic situation we 
do expect deviations from pure scale invariance induced by damping processes during the 
decoupling era (e.g. frozen amplitudes in the radiation era, neutrino free streaming or 
photon diffusion in the case of Silk damping). In numerical simulations the cutoff is intro- 
duced either by the grid or discreetness effects in gridless N-body simulations. In the latter 
case, l/k c or r in equation (32) would correspond to the mean interparticle separation (the 
particles Nyquist scale). In any case, all this suggests that scale-free power spectra are just 
convenient analytical approximations and the need for a cutoff in some perturbation theory 
applications may just tell us when this approximation has to break down. Of course, in a 
real case, the regularization may come by a smooth bending on the spectra rather than a 
sharp cutoff. 

A cutoff of this kind affects the initial conditions for the field 61, i.e. the properties of 
the density field itself. In contrast, spatial smoothing applied at late times, as in equation 
(27), describes a particular way of measuring the variance, (S 2 ). These two kinds of window 
functions should not be confused with each other. In linear perturbation theory, for a given 
choice of a window and a power spectrum, there is no difference between a field that has 
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been spatially smoothed at some initial time and then evolved, and a field which was 
evolved first and smoothed later. This is no longer true when weakly nonlinear effects are 
taken into account: smoothing and evolving do not commute. To keep the distinction, 
we reserve a capital R for the comoving smoothing radius of a filter applied at late times, 
while lower-case letters, r and k c , for cutoffs present in the initial conditions. 



3 Results and discussion 
3.1 Scale-free power spectra 

Consider a power spectrum of the form given by equation (31). In the case of no smoothing 
(R = 0, W = 1), the linear variance and the nonlinear correction term both diverge when 
the cutoff is removed by setting k c = oo. However, when the correction term and the linear 
variance are calculated for a finite k c , their ratio, I2, remains finite in the limit k c — > 00. 
By numerically integrating J 2 2 and ii 3 with W(kR) = 1, we obtain 

I 2 = 1.83 . (33) 

This value is very weakly dependent on n (contrary to the lowest order values of normalized 
cumulants). It decreases slightly with growing spectral index n, but differences in the range 
—3 < n < 1 are smaller than 1 %. Therefore we may approximate the weakly nonlinear 
(<5 2 ) as 

°lnl = ° 2 + 1-83 a 4 . (34) 

The above relationship is in excellent agreement with a recent result, obtained indepen- 
dently by Scoccimarro and Frieman (1995), who succeeded in evaluating the integrals 7 2 2 
and ii3 analytically, and obtained 

4007 

h = /22 + /13 = 2205 ~ L82 • (35) 

This expression ignores a small n-dependent term, which can affect the sum of J 2 2 and ii 3 
by no more than 1% in the entire range — 2 < n < 2 (again, in agreement with our results 
of numerical integration). 

As we have no other way of checking whether the perturbative series does converge to 
the true solution, it is important to compare perturbative results with N-body simulations. 
For the results in equations (34) and (35) this is not possible because they assume no 
smoothing (R = 0), while any experimental measurements in N-body simulations, as well 
as in galaxy surveys, necessarily involve spatial smoothing over some finite scale R ^ 0. 
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In order to account for spatial smoothing, we have calculated numerically J 2 = -^22 + ha 
for a range of spectral slopes, using a Gaussian filter (eq. [16]). Figure 1 shows the results 
of the calculation of the weakly nonlinear corrections (<5f) + 2 (S183), divided by We 
used power-law spectra with a shortwave cutoff at k = k c (eq. [31]) (therefore the value of 
a 2 G is not the one given by eq. [22] but smaller, dependent on the cutoff). The corrections 
are shown for different values of the dimensionless variable k c R. We see that in the limit 
k c — > 00, the integrals are convergent only when n = —2. However, for any finite value 
of k c R, their sum, 1 2 remains finite, and it is easy to estimate its value analytically for 
k c R ^> 1. Indeed, we can use the closed form expressions for P 2 (^) = -P22 + -P13, obtained 
by Makino et al. (1992), who studied quadratic corrections for power-law spectra with 
shortwave cutoffs. (Although their formulas are not free from errors, the mistakes we 
identified do not affect the results presented below.) One can then identify the terms, 
which dominate in the limit k c — > 00, substitute the result in equation (27), and then 
integrate term by term. For n = — 2, the appropriate expansion for P 2 is 



c 2 



(2vr) s 



55vr 2 
196 A; 



This gives 



for a Gaussian window, and 



557T 

196 

1375 
1568 



+ 0(l/k c ) 



0.882 



0.877 



(36) 



(37) 



(38) 



for a top-hat. The same approach can be used for the remaining cases, considered in 
Figure 1. Upon expanding P 2 (^) m powers of k c and neglecting the terms of order l/k c 
and higher, one has 



P2(k) 



C 2 



c 2 



(2vr) 2 
C 2 



1/2 k\n(k/k c ) + -!^-k+ 0(1/ k c ) 



315 



122 



11025 
lOvr 2 



k 2 k c + k 3 + 0(l/k e ) 

315 147 v 1 ' 



(n = -l); 
(n = 0) ; 



4973 8 



k k„ ~\~ k kp 
315 c 49 V44100 105 



ln(Jfc/A; c ) k 5 + 0(l/k c 



(39) 

(40) 

(n = l) 
(41) 



Substituting above expansions into equation (27) with a Gaussian window, and integrating 
termwise, we finally get 



I 2 (k) = 0.638 - 0.387 ln(k c R) (n = -1) ; 
h (k) = 1.710 - 0.656 k c R (n = 0) ; 



(42) 
(43) 
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J 2 (k) = -0.321 + 0.457 \n(k c R) + 0.610 k c R - 0.387 (fc c i?) 2 (ra = 1) . (44) 

For k c R > 2, all of the above asymptotic expressions for I 2 , as well as the estimate of J 2 for 
n = — 2 in eq. (37), are in excellent agreement with the results of numerical integration, 
shown in Figure 1. 

Let us now discuss the limits of applicability of the second order perturbation theory 
as a function of k c R. We must be particularly careful in the limit of large k c R, since 
(S 2 ) = a 2 [ 1 + a 2 I 2 ] , while 

|7 2 | « (0.6fc CJ R)" +1 for n > -1 and k c R^>l. (45) 

Hence, the term quadratic in variance could "overtake" the linear term, causing a break- 
down of perturbation theory. From these considerations we expect that the approximate 
limit of validity of the second order results on a scale R is given by the condition 

a(R) < \h\-^ 2 . (46) 

The upper limit on k c R is therefore 

k c R < (i?/^) ( " +3)/(n+1) (n > -1) , (47) 

where R n i is the filtering scale, for which the rms density fluctuation equals unity: a(R n i) = 
1. As we have seen, the results for n = — 2 are not sensitive to the upper cutoff k c , so the 
limit of their validity is given by the usual condition a < 1 - a standard requirement for 
perturbation theory. Unless k c R is "exponentially large" , the same limit should apply for 
the n = — 1 case, since the divergence in this case is only logarithmic. However, for n > — 1, 
the condition (46) becomes more and more restrictive with growing n. In particular, we 
expect that the range for acceptable values of a for n = 1 is narrower than the range for 
n = 0. 

To check the self consistency of our results, we will now consider the limit R — > oo 
for a fixed upper cutoff k c and a fixed nonlinear scale R n i. We should recover the linear 
perturbation theory, which is indeed the case for n > — 1 : 

For — 3 < n < — 1 the convergence to linear theory is even faster. 

In the opposite limit, when R = 0, we expect to recover the result for J 2 without 
smoothing, since W(0) = 1. Indeed, for k c R = 0, all curves in Fig. 1 do converge on 
J 2 = 1.83, in agreement with equation (34). 



(48) 
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The ultimate test for perturbative calculations is the comparison with N-body experi- 
ments. We will make such a comparison in the next section. Meanwhile, let us summarize 
the physical implications of equations (37) and (42) - (44). According to these results, 
weakly nonlinear effects accelerate the growth rate of density fluctuations when n = — 2: 
the rms density contrast grows faster than predicted by linear perturbation theory. How- 
ever, when there is more power at small scales in the initial spectrum, the nonlinear effects 
work in the opposite direction and the growth is slowed down. For n = — 1, the quadratic 
correction to linear variance is close to zero - the rate is in almost perfect agreement with 
linear theory. For n > —1, the growth rate is slower than the linear prediction. The 
quadratic correction has a negative sign and its absolute value appears to increase with 
increasing n as (k c R) n+1 , at least for the cases n = and n = 1, considered here. This 
behaviour is in qualitative agreement with the previrialization conjecture (Peebles & Groth 
1976; Davis & Peebles 1977). 

3.1.1 Comparison with N-body simulations 

To compare predictions of perturbation theory to the N-body results in the case of spectral 
indices n = — 1, n = and n = +1 we use the simulations made by David Weinberg that 
were used to check the perturbative calculation of skewness and kurtosis (Juszkiewicz et 
al. 1995, Lokas et al. 1995). All the simulations used a 200 3 force mesh and 100 3 particles 
(except for n=+l ones which had 200 3 particles). The moments of the evolved density 
field were computed for Gaussian smoothing lengths L/50, L/25 and L/12.5, L = 100 cells 
being the size of the simulation box. 

Figure 2 compares the N-body results to the perturbative calculations. Open symbols 
show the ratios of the N-body nonlinear variance to its linear counterpart calculated from 
equation (22) using the normalization of the initial power spectra (i.e. a = 1 for the 
final expansion factor a = 1 and smoothing scale L/50). Circles, triangles and squares 
correspond to the shortest, medium and largest smoothing scale respectively. The error 
bars of the results (not plotted) coming from statistical averaging over eight independent 
simulations are large enough (especially in the n = — 1 case where the points are most 
scattered) so that the results do not contradict self-similarity. 

A direct, quantitative comparison between the perturbative and N-body results for 
n > —2 is restrained by the difficulty of determining the real small scale cutoff in the 
initial power spectra of the simulations. The degree of self-similarity displayed by the N- 
body results would indicate that the cutoff scale is very small. We performed perturbative 
calculations of the weakly nonlinear corrections to the variance for the spectrum (32) with 
a Gaussian cutoff. The effect of a Gaussian cutoff in A;-space is similar to that of a sharp 
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cutoff at k = k c ~ 1/r (eq. [32]). We chose three different cutoff wavelengths r: r^y, f Ny/^ 
and r Ny /4:, where r Ny is the scale corresponding to the Nyquist frequency. The Nyquist 
scale is L/50 for n = — 1 and n = simulations and L/100 in the case of n = +1. The 
choice of a cutoff scale r obviously breaks the self-similarity of the results. We showed in 
Figure 2 only the results of perturbative calculation corresponding to the medium final 
smoothing radius L/25 (filled triangles). We find that the perturbative results match the 
N-body best when we assume the cutoff scales r = r Ny /k for n = —1, r = r Ny /2 for n = 
and r = r Ny for n = +1. These values correspond to k c R = 8 for n = — 1 and k c R = 4 
for n = and n = 1. Using equations (42), (45) and (46), we get a < 1, a < 0.6 and 
a < 0.4 as approximate limits of validity of perturbation theory for n = —1, and +1, 
respectively. These limits, based on the comparison of the linear and quadratic terms, are 
in reasonable agreement with the true limits of validity, implied by Figure 2. 

In the n = — 2 case, we performed a new simulation with 256 3 particles (about 17 
million) using a PM code (Bouchet, Adam k Pellat 1985; Moutarde et al. 1991) with 256 3 
cells. The initial conditions were imprinted using Zeldovich approximation on a "glass" -like 
particle distribution (White 1994), with a power equal to 1/25 of the shot noise level at the 
Nyquist frequency of the particle grid. This "glass" distribution was obtained by N-body 
simulation in an expanding universe with a repulsive gravitation, starting from a random 
initial distribution. After a high enough expansion (about 10 6 in our case), particles settle 
down in a quasi-equilibrium state. This state shows a very uniform particle load without 
any anisotropy or discernible order down to very small scales. This homogeneity allows 
a very accurate study of structure formation in high density region, as well as in voids. 
The variance of the evolved density field was computed on the discreet distribution of 
a 64 3 particle sub-sample convolved with a spherical top-hat of radius R. Poisson noise 
associated with discreetness effects has been removed from measured variances. 

Figure 3 compares the perturbative predictions and the N-body results for n = —2. 
Open symbols show the ratios of the nonlinear variances to the linear ones, at different 
expansion factors and at 8 scales, starting from R = 10~ 2 2 times the box size, and spaced 
by 0.2 in log. The linear variances were calculated according to equation (23) and the 
perturbative weakly nonlinear approximations to the nonlinear ones (filled symbols in 
Figure 3) using equation (38). The N-body results closely follow theoretical predictions and 
obey a self-similar evolution, apart from the largest scale measurements (with R = 10~ 8 
times the box size). These deviations are likely to be due to the effect of the missing waves 
at scales larger than the numerical box. The absence of these long waves can also explain 
the slight offset of the numerical results vs. the theory, because according to Figure 1, 
for n = —2, the correction to linear theory for variance is actually quite sensitive to long 
waves contribution. We have calculated the perturbative corrections with a cutoff at low 
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wavenumbers corresponding to the size of the box for a few points in Figure 3 and found 
that introducing this cutoff decreases both the linear variance a 2 and the weakly nonlinear 
approximation of (<5 2 ), but their ratio also is decreased. Therefore the effect of such a 
cutoff is generally to decrease the perturbative values in Figure 3 so that the theoretical 
curve moves closer to the N-body results. 

We can summarize this section as follows. All our N-body experiments support per- 
turbation theory: nonlinear evolution makes ( S 2 ) larger than the linear variance, when 
n = —2. There is an opposite effect - a slow-down of growth rate, when n > —1, and 
the magnitude of the slow-down (the value of J 2 , determined from N-body experiments) 
increases with n. The transition occurs at n = — 1, when the nonlinear correction term is 
close to zero. 



3.1.2 Comparison with the Hamilton ansatz 

Hamilton et al. (1991) proposed a general formula relating the linear two-point correlation 
function of arbitrary shape and its strongly nonlinear counterpart. While physically moti- 
vated in its limits, the overall functional shape of the relation was modeled using N-body 
simulations of scale-free initial power spectra with fl = 1. A similar formula for power 
spectra was obtained by Peacock & Dodds (1994) and recently refined to take into account 
the dependence on the spectral index of the spectrum by Mo, Jain & White (1995). We 
will use this n-dependent formula to calculate the nonlinear variance and compare the 
result to the linear prediction as was done in the case of perturbative approximation and 
N-body simulations. 

The ansatz for the relation between the linear power spectrum (A L (k ) = A.nklP{kQ)) 
and the nonlinear (evolved) one (A E (k) = A.nk 3 P{k)) is 



A L (k ) 
B{n) 



(49) 



B(n) 

where the wavevectors k and k are related by 

k = [1 + A E (k)]^ 3 k (50) 

and B(n) is a constant depending on the power spectrum index n. The values of B{n) are 

1.64, 1, 0.54 and 0.24 for n = +1, n = 0, n = — 1 and n = — 2 respectively. We use the 
fitting formula of Mo et al. (1995) 

[l + 2x 2 -0.6x 3 -l.bx 7 ' 2 + x*\ 1/2 . , 

= x (51) 

V ' \ l + 0.0037x 3 ) v ; 
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to calculate the nonlinear power spectrum A E at specified values of ho- Then, using 
equation (50) we find the value of k corresponding to each of pairs [k , A E (k )]. The list of 
points [k, (A;)] is then fitted numerically by the linear power spectrum plus polynomial 
terms of higher order. This fitted shape of the nonlinear power spectrum A E (k) is then used 
to calculate the nonlinear variance < S 2 >. For simplicity we use a filter function in the 
form of a top-hat in Fourier space (that is we cut off the integration at some wavenumber 
k c ), the effect of which should be close to the effect of a Gaussian filter with the smoothing 
scale corresponding to k c . The comparison of the results with the linear variance, a 2 , is 
shown in Figure 4. Given the numerous approximations applied in the calculations, the 
agreement between the results thus obtained with those of perturbative calculations and 
N-body simulations (Figures 2 and 3) is rather impressive. 



3.2 The Peacock-Dodds spectrum 

Scale-invariant spectra can be regarded as idealizations of "realistic" spectra, which arise 
either from a mixture of model assumptions and physics of the decoupling era, or from 
attempts to parametrize the observational data. Such spectra usually have slopes depen- 
dent on the wavenumber, n(k) = d (logP) / d (logfc), therefore the size and sign of weakly 
nonlinear corrections to P{k) can be also expected to vary with k. As an example of such 
a spectrum we consider the following 

A(k) - {k/k ° r (52) 

where A(k) = i-n:k 3 P(k) /(2tt) 3 . The parameters given by Peacock & Dodds (1994), 

k = 0.29 ± 0.01 h Mpc" 1 

k c = 0.039 ± 0.002 h Mpc" 1 (53) 
a = 1.50 ±0.03 
P = 4.0 ±0.5, 

were fitted to best match the data obtained from eight independent galaxy surveys. The 
reconstruction of this linear spectrum involved accounting for nonlinear evolution as well 
as redshift space distortions. If we accept the most probable value a — (3 = — 2. 5 we may 
rewrite the spectrum in the following way 

Ck 

P{k) = ^L— (54) 

w 1 + (k/k c ^ n 



with C = 2% 2 k x ' 5 k c 2 5 and n = 2.5. The spectrum of this kind was first considered by 
Peacock (1991) who quotes k c in the range [0.015, 0.025] h Mpc" 1 and n = 2.4. 
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Following equation (19), the linear variance of the density contrast averaged over a 
sphere of radius R is 

2 9 71 " [°° d& 2 

° R = Wj¥r 4 'o i + (^) 2 - 5 3/2 ( ) ' (55) 

which is easily integrable numerically. Assuming the smoothing scale R = 8 h~ l Mpc we 
obtain 

a s = 0.69. (56) 

We calculated the nonlinear correction, using equations (27), (28) and (29) and a convenient 
approximation of a top-hat window function, 

W TH {kR) « exp (-k 2 R 2 /9) , (57) 

which is an accurate representation of a real top-hat (eq. [17]) to within few percent in 
the entire range of integration. The convergence of the integrals is now ensured because 
for k — > oo the power spectrum (54) approaches C/k 1 ' 5 . The calculation must however be 
done for every considered smoothing scale independently as the spectrum is not scale-free. 
At R = 8 h- 1 Mpc we find 



7 2 (8/ l - 1 Mpc) = ^ + 2 W 8 = 0.15 



(58) 

°8 

and the value of a 8 taking into account the weakly nonlinear corrections is 

<7 8jU)nI = 0.72 (59) 

which is less than 4 % above the linear value. This result is consistent with what we have 
shown above for the scale-free power spectra. At the scale of 8 h~ l Mpc the spectrum 
has index close to —1 or slightly below this value. For such spectral indices we expect the 
weakly nonlinear variance to be equal to or slightly above the linear value. 

One may ask how sensitive these results are to uncertainties in fitting parameters, given 
in equations (53). If we take the claimed error bars seriously, the resulting uncertainty in 
I 2 in equation (58) is rather small. Varying the least accurate parameter, j3, from 3.5 to 4.5 
gives 7 2 = 0.15 ± 0.02. Varying the other fitting parameters in the allowed ranges leads to 
negligible changes in I 2 . Of course, we must keep in mind that the true uncertainties in the 
shape of P(k) may be quite different. The parameters in equation (53) were not obtained 
from a homogeneous data set. Instead, the process of data reduction required a number 
of approximations, including a procedure applied to unify the data coming from different 
surveys and the method to account for the nonlinear evolution which is not expected to 
work equally well for different parts of the power spectrum. 
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3.3 Concluding remarks 



We have computed the first nonlinear corrections to the variance of a density field under- 
going gravitational instability. The results are confirmed by numerical simulations, and 
agree with the ansatz of Hamilton et al. (1991), as modified by Mo et al. (1995). There 
is indeed a previrialization effect, as conjectured by Peebles & Groth (1976) and Davis & 
Peebles (1977) and Peebles (1990). The effect depends on the initial power spectrum con- 
sidered. For n < — 1, the fluctuations grow faster than in linear theory, while for n > — 1, 
we see "previrialization" , or slow-down. 

The magnitude of this effect grows when the relative amount of small-scale power is 
increased by increasing n. The effect is small at n = — 1, corresponding to the transition 
between the slow-down and speed-up of the collapse. Our results can be used to explain 
the possible cause of differences of opinions about the size of the previrialization effect. 
Indeed, Evrard & Crone (1992) assumed n = — 1 in their simulations, and saw no effect, 
while Peebles (1990) reached an opposite conclusion, using simulations with n = 0. 

For scale-free power spectra the difference between the linear and weakly nonlinear 
approximation for the variance can be as high as 100% as in the case of the spectral index 
n = —2. For the realistic power spectrum of a class considered by Peacock & Dodds (1994) 
however, we have found that the correction induced by weakly nonlinear effects on a 8 is 
very weak. This is purely by chance: the effective index of the realistic spectrum at the 
scale R = Mpc happens to be close to n = — 1. For such an index, as we have seen, 
the nonlinear correction is close to zero. 

Our results are in agreement with previous studies of nonlinear interaction between 
perturbations at different scales. Using the fluid model for the evolution of structure 
Peebles (1987) found that for the initial power spectrum of index n = (and some small 
scale cutoff) the smoothed standard deviation of the evolved field is smaller than its linear 
extrapolation. Weinberg & Cole (1992) performed a series of N-body simulations with 
Gaussian initial conditions and scale-free initial power spectra normalized so that the 
evolved cr 8 = 1. For the power spectrum index n = — 1 they found that a 8 grows at almost 
exactly the rate predicted by linear theory, while for n = the required linear cr 8 was 
larger than the evolved value and for n = — 2 smaller. Jain & Bertschinger (1994) found 
that for CDM spectrum normalized so that linear cr 8 = 1 the second-order effects increase 
a 8 by 10%. 

Why is the n = — 1 index so special? At least a partial explanation may be provided 
by recalling the properties of the peculiar gravitational acceleration, g , known from linear 
perturbation theory. For a density field, smoothed on scale R, g oc i? _ ( n+1 )/ 2 (Peebles & 
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Groth 1976; Vittorio & Juszkiewicz 1987). If n < —1, g diverges at large R: the peculiar 
acceleration field has an infinite coherence length (the same is true for the velocity field). 
The gravitational acceleration moves the fluid more or less uniformly, generating bulk 
flows rather than shearing motions. Such an acceleration field will move a collapsing mass 
concentration as a whole, without affecting its internal dynamics. Therefore, its collapse 
will be similar to that of an isolated spherical clump (in which 6 grows faster than in linear 
theory). In the opposite regime, when n > — 1, the dominant sources of acceleration are 
local, small-scale inhomogeneities. Hence, we can expect a slow-down of the growth rate: 
tidal effects will tend to generate nonradial motions and resist gravitational collapse. The 
above picture, based on linear theory arguments, can be made more rigorous by involving 
higher terms in the perturbative expansion and by studying the evolution of high density 
peaks rather than the evolution of regions, selected at random. This was recently done by 
Bernardeau (1994). The transition near n = — 1 was also studied in N-body simulations 
by Klypin & Melott (1992) and Melott & Shandarin (1993), who concentrated on the 
properties of the peculiar velocity field. 
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Figure captions 



Figure 1 Weakly nonlinear corrections to the variance (<5f) + 2 (6163} , divided by cr^, 
(eq. [30]) for a Gaussian filter and different scale-free power spectra. Symbols show 
the results of numerical integration up to different values of k c R, where k c is the 
cutoff wavenumber in the initial spectrum and R is the scale of the final smoothing. 
Only for n = — 2 the integration converges to 0.86, which is slightly different from 
the value obtained analytically (eq. [37]). In the limit of k c R — > all curves converge 
to values corresponding to the case of no smoothing (eq. [33]). 

Figure 2 The ratio of the nonlinear to linear variance, (<5 2 ) / a 2 , as a function of a 2 for 
different initial power spectra. Open symbols correspond to the results of N-body 
simulations with Gaussian smoothing lengths of L/12.5 (squares), L/25 (triangles) 
and L/50 (circles), where L is the size of the simulation box. Filled triangles show 
the results of perturbative calculations with final smoothing radius L/25 and the 
small scale cutoff r^y/A for n = —1, r^y/2 for n = and r^ y for n = +1. 

Figure 3 The comparison of the perturbative vs. N-body results for the power spectrum 
of index n = — 2 and top-hat smoothing. Open symbols show the ratio of the 
nonlinear to linear variance, measured in N-body experiment at different expansion 
factors and eight final smoothing scales starting from R = 10~ 2-2 times the size of 
the simulation box and spaced by 0.2 in log. The corresponding perturbative results 
(filled symbols) were calculated assuming the normalization of the simulated initial 
spectrum and using equation (38). 

Figure 4 The ratio of the nonlinear to linear variance, as a function of a 2 for different 
initial scale-free power spectra. The nonlinear variance was obtained from the non- 
linear power spectrum given by an ansatz, similar to the one, proposed by Hamilton 
et al. (eqs. [49] - [51]). Both linear and nonlinear values were calculated using a 
top-hat filter in Fourier space. 
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